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ABSTRACT 


Results  on  the  statistical  analyses  of  series  of  events  published 
subsequent  to  the  monograph  by  Cox  and  Lewis  on  this  subject  are  surveyed, 
Special  emphasis  is  given  to  tests  for  renewal  processes,  a  considerable 
amount  being  now  known  about  the  distributions  of  some  of  the  test  sta- 
tistics involved,  and  to  testing  the  functional  form  of  a  trend  in  a 
nonhomogeneous  Poisson  process,  as  well  as  the  point  process  model  itself 
A  survey  of  work  in  special  processes  such  as  cluster  processes  and 
doubly  stochastic  Poisson  processes  is  also  given. 

Prepared  by: 


1.     INTRODUCTION 

The  object  of  this  paper  is  to  survey  recent  results  in  the  statisti- 
cal  analysis  of  univariate  point  processes  (series  of  events).      The  survey 
is  personal   reflecting  my  own  present  interests,   and  is  not  comprehensive, 
For  convenience,    I  have  taken  "recent"  to  mean     anything  published  since 
Cox  and  Lewis  (1966).     I  have  also  mentioned  some  topics  and  results 
which  were  omitted,   for  reasons  of  emphasis  or  ignorance,   from  that 
monograph.      Finally,   I  point  out  areas  where  further  work  is  required. 

A  survey  of  work  in  this  field  is  not  so  difficult  as  a  survey  of, 
say,   the  theory  of  point  processes  because  the  advances  in  the  statisti- 
cal analysis  of  point  processes  are,   by  comparison  with  the  theory,    few. 
This  may  reflect  the  relative  difficulties  of  these  areas,   but  that  thought 
may  be  a  personal  bias. 

Some  of  the  shortcomings  of  the  Cox  and  Lewis  monograph  were, 
in  hindsight,   the  following: 

(i)  Not  enough  consideration  of  grouped  data.     In  some  cases  where 

point  events  occur,    such  as  epidemiology,    recording  limitations  or  the 
volume  of  data  force  one  to  work  with  numbers  of  events  in  fixed  inter- 
vals whose  width  may  or  may  not  be  controllable  in  advance.     I  have 
touched  on  this  problem  recently  (Lewis,    1970);     separate  spectral 


analyses  of  the  intervals  and  of  the  counting  process  (Bartlett,    1963) 
are,   for  example,   not  possible.     One  has  to  use  a  spectral  analysis  of 
the  grouped  counts    (i.  e.  ,   number  of  events  in  a  fixed  time  interval), 
which  is  very  nonstandard  in  its  distribution  theory,   but  in  which  well 
known  problems  of  aliasing  arise.     Results  of  Cox  (1970)  may  be  useful 
here,    and  Cox  (1955)  has  considered  some  other  aspects  of  the  analysis 
of  grouped  events. 

Interesting  examples  of  this  type  of  problem  are  found  in  physics 
and  optics.     For  example   (Helstrom,    1964;  Karp  and  Clark,    1970) 
photon  or  other  particle  emissions  are  known  from  physical  consider- 
ations to  be  generated  by  a  doubly  stochastic  Poisson  process  and  it  is 
required  to  determine  parameter  values  of  the  driving  process  from 
counts  of  the  number  of  photons  emitted  in  successive  periods.      Here 
it  is  prohibitively  costly  to  record  exact  times  of  occurrence.     However, 
the  recording  interval  can  be  determined  in  advance  by  the  experimenter. 
Problems  of  grouped  Poisson  counts  (McNeill,    1971)  are  also  common. 

(ii)  Very  little  emphasis  was  given  to  sequential  methods.     While  this 

was  to  deliberate  to  save  space,   it  is  also  true  that  most  data  I  have  come 
across  is  presented  for  analysis  in  toto.      This  may  change  as  better    re- 
cording methods  are  introduced  and,   hopefully,   as  statisticians  are  called 
in  before  the  fact.      There  is  still  a  problem  in  that  simple  sequential 


methods  are  known  only  for  Poisson  processes    (homogeneous  or  non- 
homogeneous)    and  also  that  rather  unsmooth  inhomogeneities  occur  in 
practice  which  make  the  application  of  formal  sequential  methods  based 
on  very  definite  assumptions  quite  hazardous. 

Less  formal  sequential  methods  are  useful;    in  particular,    an 
analysis  of  the  data  in  successive  sections  is  very  useful,    both  to  cut 
down  on  computation  time  in,    say,    spectral  analyses     (Lewis,    1970) 
and  to  examine  the  time  evolution    of  the  process. 

As  an  example,    consider  a  series  of  arrivals  at  an  intensive 
care  unit  in  a  hospital.      This  data  will  be  used  for  illustration  through- 
out the  paper;  it  was  supplied  by  Dr.    A,    Barr,   of  the  Oxford  Region- 
al  Hospital  Board,    England.  The  first  section,    consisting  of 

n=  251    arrivals  in    t     =  409    days    (4  February,    1963  -  18  March,    1964), 
was  analyzed  in  Chapter  1  of  Cox  and  Lewis    (196  6).      Later  on,    the 
arrivals  up  to    6  February,    1968    were  received.      Three  subsequent 
sections  of  length    t     =  409    days  were  taken  from  these  later  arrivals 
for  comparison  and  their  statistics,    as  well  as  those  for  the  total 
record,    are  shown  in  Table  1. 


*    These  times -of -arrival  were  exceptionally  well  recorded.     Of  the 
146  8  arrivals  in  the    142  0-day  period  from    4  February,    1963    to 
6  February,    196  8,     only  one  time-of-day  was  not  recorded.     Nine  tied 
arrivals  occurred.     Generally,    recording  times  seemed  to  be  at  the 
five -minute  intervals  of  the  hour,    although  other  times  do  occur.      The 
data  to    18  March,    1964    is  given  in  Cox  and  Lewis  (1966,   p.  255);     the 
rest  in  Lewis     (1971). 
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Table  1.  Arrivals  at  an  intensive  care  unit  (no  ties) 

Period  1  Period  2             Period  3            Period  4           Total 

4  Feb  '63-  19  Mar  «64-       3  May  '65-       16  June  '66-    4  Feb  '63- 

18  Mar  '64  2  May  '65           15  June  '66       29  July  '67       6  Feb  «68 

n                   251  350 

fc0                 409  409 

U(days)       218.47  218.53 

U                       1.874  2.222 

msN/t.          0.  614  0.  856 


n 
The  statistic     U  =    2    t./n,     where  the    t.     are  the  timet  to  events  and 

i«i  *  * 

n    the  number  of  events  in  the  period  of  observation  of  length    t   ,     is  the 
centroid  of  the    t. 's    and  is  used  to  test    (Cox  and  Lewis,   1966,   p.  47) 
for  the  trend  parameter    p     in  a  nonhomogeneous  Poisson  process  with 
rate 


X  (t)  «  exp(er  +  pt)*  X.    exp(pt).  (1.  1) 

In  testing  for     p  *  0    against    p  ^  0,     or    is  a  nuisance  parameter  with 
sufficient  statistic    n    and  the  test  is  based  on  the  distribution  of    U, 
given    n.     The  normalized  statistic    U    in  line  three  of  Table  1      should 
be  distributed  approximately  as  a    N(0,  1)    variable.      There  was  fairly 
strong  evidence  of  a  monotone  trend  at  the  end  of  the  first  observation 
period    (  P*  7  percent  level),     but  later  on  in  the  series,    as  in  Section  3, 
there  is  a  definite  decreasing  trend.     However,    the  total  arrival  process 
gives  fairly  strong  evidence  of  an  increase  in    X  (t),     the  value    U  ■  2.  874 


being  significant  at  about  a  2  percent  level.    For  Period  1   and  Period  2  combined, 
U   was  found  to  be  4.492  which  is  significant  at  a  level  much  smaller  than  1  percent. 
An  actual  sequential  test  of    p  =  0    against    p   ^  0    rejects  the  null  hypo- 
thesis at  a    1  percent    level  after     550    days.     The  sectioned  analysis  of 
four  periods  suggests  a  long  cycle  or  quadratic  term  in  the  exponential 
trend    (1.  1).      The  nonhomogeneity  is  confirmed  very  strongly  by  a  dis- 
persion test    (Cox  and  Lewis,    1966,   p.  232)    applied  to  the  numbers  of 
events  in  the  four  sections.      This  has  value    26.  74;    its  distribution  is 
that  of  a  chi-square  variable  with    3    degrees  of  freedom  which  has  a 
0.  9999    point  equal  to    21.  11. 

A  plot  of  the  times -to-events    t.     against  serial  number    i    is 
given  in  Figure  1,    and  a  plot  of  the  average  values  of  successive  sets 
of  twenty  intervals  between  arrivals  is  shown  in  Figure  2.     Again  the 
long  cycle  or  quadratic  trend  is  quite  evident  graphically.      We  return 
to  this  example  later  in  the  paper. 

(iii)         Sectioning  brings  up  in  some  ways  the  analysis  of  replicated 
point  processes,    which  again  was  not  considered  in  detail  by  Cox 
and  Lewis  (1966).      This  replication  can  occur  quite  commonly  in  exper- 
imental situations,   for  instance  in  neurophysiology  where  experiments 
can  be  repeated  many  times.     Here  the  signals  are  trains  of  very  narrow 
spikes  of  apparently  fixed  height,    so  it  is  appropriate  to  analyze  the 
times  of  occurrence  of  the  spikes  as  a  point  process.      Observation  over 
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Figure  1.     Arrival  of  patients  at  an  intensive  care  unit. 

Complete  record,    4  February  1963  to  8  February  1968.      Number  of  arrivals    (n)  1458    in 

periodof    tn  «   1829  days.      Arrival  number  vs.    time  to  arrival,      m  =  n/t„  =  0.  7972. 
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Figure  2.      Arrival  of  patients  at  an  intensive  care  unit. 
Averages  of  successive  groups  of  20  inter-arrival  times.     Number  of  arrivals    (n)  1458 
in  periodof    tn  «=  1829  days.     Overall  average  is  1.  254  days.     Homogeneity  of  variance 
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18.  000  (x271:  mean  71,  v  =  '11.  92) 


too  long  a  period  can  be  deceiving,    because  of  physiological  deterioration, 
and  replication  is  sometimes  preferable.      Of  course,    care  must  still  be 
taken  that  physiological  or  experimental  conditions  have  not  changed. 

Comparison  of  rates  and  trends,   mainly  in  Poisson  processes, 
was  discussed  in  Cox  and  Lewis    (1966).     See  also  Qureishi    (1964)    for 
the  comparison  of  rates  in  two  Weibull   processes.     General  problems 
of  multiple  point  processes    (multivariate  point  processes)    are  discussed 
in  Cox  and  Lewis    (197  2)    and  Perkel,   Ger stein  and  Moore      (1967b),    but 
are  beyond  the  scope  of  this  paper.      Lewis     (197  0)    discussed  estimates 
of  the  spectrum  of  counts  from  sectioned  data,    as  well  as  pooling  and 
comparison  of  the  spectra. 


(iv)  Although  trend  analysis  was  discussed  in  Chapter  3  of  Cox 

and  Lewis  (1966),    I  don't  believe  it  received  enough  emphasis. 
Also,   problems  such  as  the  analysis  of  cyclic  trends  were  barely 
touched  upon.      These  and  other  problems  in  trend  analysis  are 
discussed  in  more  detail  later  in  the  paper. 


(v)  Finally,   perhaps  more  emphasis  could  have  been  given  to  graph- 

ical methods.      These  were  discussed  in  Chapter  1,    but  no  reference  was 
made  to  probability  plotting,    for  example,    in  later  chapters.      See,   for 
example,    Wilk  and  Gnanadesikan    (196  8). 


The  field  where  moit  use  has  been  made  of  techniques  for  the 
analysis  of  point  processes  is  in  neurophysiological  work  on  the  signals 
occurring  on  nerve  fibers.      A  summary  of  techniques  for  this  type  of 
analysis  given  in  Perkel,   Gerstein  and  Moore    (1967a,    1967b)    largely 
parallels  Cox  and  Lewis     (1966),     with  more  information  on  the  special 
problems  of  neurophysiology.     For  later  work  and  applications,    see, 
for  example,    Moore,    Segundo,    Perkel  and  Levitan    (197  0)    and  the 
references  given  in  that  paper.     It  would  be  impossible  to  try  to  sum- 
marize all  of  this  work  here;     some  of  it  will  be  touched  on  later,    but 
there  is  virtually  no  statistical  methodology  given  in  them  which  is  not 
given  in  Cox  and  Lewis     (1966). 

Another  interesting  field  of  application  is  to  the  study  of  the 
occurrence  of  earthquakes.      While  strictly  not  a  univariate  point  pro- 
cess,   an  approximation  to  the  earthquake  process  as  a  univariate  point 
process  yields  useful  insights.      For  such  analyses,    see  Vere-Jones, 
Turnovsky  and  Eiby    (1964),     Vere-Jones  and  Davies     (1966),     Vere-Jones 
(1970),    and  Shlien  and  Toksoz  (1970).     The  discussion  of  Vere-Jones  (1970) 
contains  extensive  comments  on  the  problem  of  analyzing  earthquake  oc- 
currence data. 

Computational  problems  in  the  spectral  analysis  of  point  processes 
have  been  discussed  by  Lewis     (1970).     By  spectral  analysis  here  and  in 
the  rest  of  the  paper  we  mean  Bartlett's  spectral  analysis  of  the  counting 
process    N(t)    of  a  point  process    (Bartlett,    1963;     Cox  and  Lewis,    1966, 
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Chapter  5).      The  spectral  analysis  of  the  intervals  between  events  in  a 
point  process  has  been  greatly  facilitated  by  the  availability  of  the  fast 
Fourier  transform  algorithm;     see  Cooley  and  Tukey  (196  5)  and  Cooley, 
Lewis  and  Welch    (197  0)    for  details. 

Most  of  the  statistical  methodology  for  the  analysis  of  univariate 
series  of  events  given  in  Cox  and  Lewis    (1966)    has  been  implemented  in 
a  computer  program  called    SASE  IV.     For  details,    see  Lewis    (1966) 
and  Lewis,    Katcher  and  Weiss    (1969).     SASE  IV    is  a  large    FORTRAN 
program  with  graphical  output  which  is  available  from  the    IBM  Program 
Information  Department,    40  Saw  Mill    River  Road,    Hawthorne,    New  York 
10532,     as  Program  No.    360  L  130001. 

In  subsequent  sections  of  this  paper  we  discuss  first  two  central 
problems  in  the  statistical  analysis  of  point  processes,    namely  tests 
for  renewal  processes  and  tests  for  Poisson  processes.      Then  techniques 
for  specific  non- renewal  processes,  such  as  doubly  stochastic  and  cluster- 
ing Poisson  processes  are  discussed.      The  inability  to  write  down  a 
likelihood  function  makes  the  formal  analysis  of  non- renewal  point  pro- 
cesses difficult;    it  is  only  for  the  important   doubly  stochastic  Poisson 
process  that  techniques  are  beginning  to  appear. 

The  final      sections  deal  with  trend  analyses,   mostly  of  non- 
homogeneous  Poisson  processes;    first  we  consider  the  case  of  monotone 
trends,    then  the  case  of  cyclic  trends  and  their  relationship  to  the 
spectral  analysis  of  point  processes.     Following  this,    some  general 
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problems  of  trend  analysis  are  considered;     these  include  tests  for  particular* 
rate  functions  and  the  nonhomogeneous  Poisson  process  model  per  se,    and 
the  definition  of  "residuals"  in  the  analysis  of  nonhomogeneous  Poisson 
models. 
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2.  TESTS  FOR  RENEWAL  PROCESSES 

2.  1  Markov  interval  processes  and  serial  correlation  coefficients. 

A  natural  extension  of  the  renewal  process  model  is  to  processes 
with  first  order  dependence  of  intervals     (Wold,    1948;     Cox,    1955).      The 
naturalness  may  be  only  mathematical    since  point  processes 
of  this  type  have  not  been  commonly  found  in  applications,   although 
they  have  recently  been  postulated  in  neurophysiological  contexts 
(Lampard,      1968;     Walle,    et  al,      1969).     In  the  process  of  Walloe 
et  al    (1969),     the  first  order  Markov  interval  property  depends 
on  the  input  to  a  neuron  being  a  Poisson  process.     However,    since  the 
input  is  the  superposition  of  an  unknown  number  of  fairly  regular  neuron 
signals,    the  hypothesis  is  tenuous.     A  drawback  in  the  use  of  the  first 
order  Markov  interval  process  as  an  approximate  model  for  serial 
dependence  in  series  of  events  is  the  difficulty  of  obtaining  analytical 
results  on,    for  example,    the  spectrum  of  counts  or  the  variance-time  curve, 
This  difficulty  is  closely  related  to  the  fact  that  there  are  no  regenera- 
tion points  in  the  process. 

Another  problem  with  the  model  is  the  dearth  (until  recently) 
of  useful  bivariate  distribution  models.  Cox  (1955)  used  a  bivariate 
exponential  of  the  form 


f         (x;   x.    =  y)  =  >^  (y)  exp{-\  (y)x>,  (2.  1) 

Xi+1  ' 
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where 


My)*^0(i  +  \x)>  o, 


which  had  the  drawback  of  nonlinear  restrictions  on  the  parameter. 
Another  widely  used  model  for  bivariate  gamma  distributions  which 
arises  quite  naturally  has  been  discussed  by  Moran    (1967a),    Vere-Jones 
(1967),    Lampard  (196  8),    and  Gaver  (197  0).     (See  also  Griffiths,    1969, 
for  general  properties  of  bivariate  gamma  distributions.  )      These  bi- 
variate distributions  all  have  positive  correlations.     Bivariate  exponen- 
tials with  negative  correlation  have  been  derived  by  Gaver  (1972)    and 
bivariate,    negatively  correlated  intervals  have  been  observed  in  a 
neurophysiological  context  by  Walloe  et  al  (1969). 

For  »tatistical  analysis,    one  of  the  most  important  theoretical 
results  on  bivariate  exponential  distributions  is  that  of  Moran  (1967a), 
who  showed  from  more  general  results  that  the  serial  correlation 
coefficient    (or  order    1)    for  bivariate  exponential  distributions  is 
always  greater  than     -0.6649; 


E{(X.   -  E(X))(X         -  E(X)} 

p« S 1+i *   1  -    i*'-  -0.6649.        (2.2) 

1  var  (X.)  6 


Estimates  of  the  serial  correlation  and  testi  of  the  renewal 
hypothesis  against  a  first  order  Markov  interval  procesi  are  based  on 
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the  sample  serial  correlation  coefficient     p  ,     defined  as  follows. 

For  simplicity,    assume  observation  starts  with  an  event  and 

ends  with  an  event,    there  being    n    observed  intervals  between  events 

x.,  x_,  .  .  .  ,  x      with  mean    x  *  n      Ex..      Let    z.  =  x.   -  x.      Then 
I      2  n  111 

n-1 

2     z.   z.    n 
.    ,      li+l 
n       1=1 

"I*    nTT        n        I  '  <2'3> 

Before  discussing  theoretical  results  on  tests  for  renewal  pro- 
cesses  based  on    p  ,     consider  its  distribution  for  renewal  processes, 
a  great  deal  about  which  is  now  known. 

The  expected  value  of    p      under  the  renewal  hypothesis,   for  any 
distribution  of  the  intervals    x. ,     is  (Moran,    1967b) 


E(?1)  -     r-lj,  (2.4) 


and  while  its  variance  is  known  to  be  asymptotically    n         if    p    =  0,   the 
exact  variance  depends  on  the  distributions  of  the    x.    (Moran,    1967b): 

2  Sz 

.  w(r,  .    »=£+L    _     »«         K{_±_  >  *     «£*  (2.5) 

1         „2(n-l)  n(n-"  (S2f)2  n(n-1) 


Moran  (1967b)    obtained  an  approximation  by  replacing  the  expectation 
by  the  ratio  of  the  expectations  of  the  numerator  and  denominator.      This 
result  is  exact  for  normally  distributed    x. 's,     giving 
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(n-2)2 
var(p)  =     v         '  ■     .  (2.6) 

n    (n-1) 


For  exponentially  distributed   x    ,      the  Moran  approximation  gives 

i 


n  n  n 


Moran  (19*7  0)     compared  his  approximation  to  sampling  results  obtained 
by  Cox  (1967)     and  a  brief  summary  of  their  results  on  the  variances 
follows. 


(a)         The  variance  of    p      tends  to  be  smaller  than  that  for  the  normal 
case  for  positive  random  variables  with  long  tails. 


(b)  Moran1  s  approximation  works  reasonably  well  with  gamma  dis 
tributions  with  shape  parameter    k  ^  0,     n^50,     or    ki8,     n  i  10. 
Outside  these  limits  the  approximation  underestimates  the  variance 

(k  =  0    is  the  exponential  distribution). 

(c)  For  the  exponential  distribution,    a  partial  reconstruction  of 
Table  1     from  Moran  (197  0)     using  sampling  results  from  Goodman 
and  Lewis     (1972)  follows. 
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________ _     A  __ 

Observed 

Sample 

variance 

Moran's 

Difference/ 

size    n 

(simulation) 

approximation 

Difference 

Observed 

20 

0.  0419 

0.  037  1 

0.  0048 

.  115 

50 

0.  0182 

0.  0176 

0. 0010 

.  055 

100 

0.  00945 

0.  00935 

0.  00010 

.  011 

450 

0.  00219 

0.  00221 

0.  00002  0 

.  009 

(d)  For  very  skewed  distributions,    a  better  approximation  to    var(p  ) 

is  needed.      This  is  apparent  from  Table  3     for  random  samples  where 

the    x.     have  a  Weibull  distribution  with  shape  parameter   — ■  ( —  Wiebull) 

A  random  variable  with  this  Weibull  distribution  is  the  square  of  a  random 

variable  with  a  unit  exponential  distribution. 

Table  3.     Variance  of    p       in    1/2  Weibull  random  samples 

Observed 

Sample  variance  Moran's  Difference/ 

size    n  (simulation)  approximation         Difference        Observed 


50 

0.  01410 

0.  00058 

0.  01352 

.96 

100 

0.  007  82 

0.  00054 

0.  00728 

.  93 

450 

0.  00201 

0.  00047 

0. 00154 

.  77 

Returning  to  the  distribution  of    p.,     Moran  (197  0b)     showed  that 

l/2~ 
n        p      converges  in  distribution  to  a  unit  normal  variate  if  the  first 

four  moments  of  the    x.     exist.      Cox  (1966)     examined  the  distribution 
by  synthetic  sampling  for  the  case  where  the    x,     had  a  normal  distri- 
bution,  and  also  distributions  of  the  form 
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r(k+i) 


1  -x     k  .         _, 

e       x   ,  (x=   0),  (2.  8) 


with    k  =  24,    8,    1,    0,    -  1/2,     a  rectangular  distribution,    a  double 
exponential  distribution  and  a  Cauchy  distribution.      The  first  two 
moments  and  coefficients  of  skewness  and  kurtosis  were  examined. 
For  the  positive  random  variables,    the  distributions  were  generally 
positively  skewed  with  convergence  to  normality  apparently  fairly 
rapid. 

These  results  have  been  extended  for  the  case  of    x.     exponen- 
tially distributed  in  large  scale  simulations  reported  in  Goodman  and 
Lewis  (1972). 

Figure  3  shows  a  computer  plot  of  the  estimated  mean  (*), 
standard  deviation  (x)    and  coefficients  of  skewness  (+)  and  kurtosis 
(A)    of    p      as  a  function  of    n    (RHO(l,  N))    for  n  -  11(1)  120,    130,    140, 
150.      The  simulations  involved  at  least    3,  000,  000    replications  for 
each    n,      so  that  the  sampling  variances  of  the  estimates  are  small. 
The  curves  have  not  been  smoothed.      Note  that  the  skewness  is  positive, 
with  a  maximum  value  of  about    0.  32    at    n  =  35,    after  which  it  starts 
back    toward  its  asymptotic  value  of  zero.      The  kurtosis  is  small, 
going  from  positive  to  negative  at  about    n  *  35    and  then  going  back 
toward  its  asymptotic  value  of  zero  very  slowly. 

The  departure  from  normality    (positive  skewness)    is  seen 
much  more  clearly  in  Figure  4,    where  the  computer  plot  gives  esti- 
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mated  normalized  values  of,    from  top  to  bottom,    the  quantiles  at 
levels     0.001,    0.002,    0.005,    0.010,    0.020,    0.025,    0.050,    0.100, 
the  estimated  mean  plus     (n-1)     ,     and  the  quantiles  at 
levels    0.900,    0.950,    0.975,    0.980,    0.990,    0.995,    0.998,    0.999 
These  should  converge  to  the  corresponding  quantiles 
of  the  unit  normal  distribution. 

Note  two  things:  the  departures  from  normality  are  relatively- 
small  for  the  inner  quantiles  and  convergence  to  the  normal  quantiles 
is  very  slow. 

More  detailed  results  from  which  the  departure  from  normality 
and  the  slow  rate  of  convergence  can  be  assessed  are  shown  in  Table  4. 
The  quantiles  and  moments  of    p      (called    RHO(l,  N)    in  the  computer) 
for    n  =  450    are  shown  in  the  rows  marked  exponential.      Of  these 
rows,    the  first  row  is  the  actual  estimated  value  of  the  moment  or 
quantile,    the  bracketed  quantities  in  the  next  row  are  the  estimated 
sampling  standard  deviations  of  the  estimates     (5  degrees  of  freedom), 
and  the  third  row  is  the  quantile  minus    u,     all  divided  by    o .      This 
row  illustrates  the  departure  from  normality,     x        _       being    2.  039 
instead  of    1.  96  0. 

Much  more  serious  departures  occur  for  the  quantiles  of    p 
from  random  samples  with    1/2  Weibull  intervals.      These  are  given 
in  the  rows  marked    "1/2  Weibull"    in  Table  4.     Thus,    not  only  is 
Moran's  approximation  poor,    as  seen  above,    but  the  normal  approxi- 
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mation  is  considerably  off.      The  simulation  results  for  the    1/2  Weibull 
case  are  discussed  more  fully  in  Goodman  and  Lewis  (1972). 

Moran  (1967b)   showed  that  tests  for  independence  against  first  order 
serially  correlated  intervals  with  exponentially  distributed  marginals 
based  on    p      are  asymptotically  most  powerful.     He  also  conjectured 
that  this  is  true  for  Gamma  distributed  intervals.     More  formal  results 
have  been  reported  by  Yang  (197  0). 

Moran  (1967b)  also  proposed  two  modifications  of  the  serial  corre- 
lation coefficient.     One  is  obtained  (see  also  Cox,    1955)  by  replacing  the 
estimate  of  the  variance  in  the  denominator  of  (2.  3)  by  the  square  of  the 
estimate  of  the  mean  (the  mean  and  the  standard  deviation  are  equal  for  the 
exponential  distribution).      The  other  modification  is  due  to  Ogawara. 
Moran  (1967b,    197  0)  gives  the  first  two  moments  of  these  test  statistics 
and  shows  that  asymptotically  they  involve  no  loss  of  efficiency  against 
the  first-order  Markov  interval  process.     However,    Moran1  s  conjec- 
ture that  the  distribution  of  Ogawara  statistic  converges  rapidly  to  the 
normal  is  not  correct    (Goodman  and  Lewis,    1972);    in  fact,    its  distri- 
bution  is  very  similar  to  that  of    p      in  the  null  exponential  case. 

An  interesting  application  of  serial  correlation  statistics  in 
examining  neurophysiological  models  is  given  by  Walloe  et  al  (1969). 
There  is  very  evident  negative  correlation  between  successive  inter- 
vals; given  the  known  structure  of  the  neuronal  process,    first  order 
Markov  dependence  would  have  been  expected  if  the  input  to  the 
neuron  had  been  a  Poisson  process,    and  this  was  tested  by  checking 
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~   2 
if    (p  )      ~  p    .      Unfortunately,      (personal  communication)     p      was  esti- 
mated by  averaging  together  estimates  obtained  from  successive  sections 

of  length  ?0,   and  the  bias  (2.  4)  accounts  for  a  large  part  of  the  observation 

i**  \  2 

that  for  several  experiments    (p    )      was  systematically  significantly  larger 

than    p    .     A  jackknifed  estimate  (Quenouille,    1956)  could  have  been  used, 
or  a  correction  for  bias  introduced.      There  probably  still  is  real  higher 
order  dependence  in  this  data  accounted  for  by  the  fact  that  the  input  to  the 
neuron  is  a  superposition  of  a  finite  number  of  inputs  and  therefore  not  quite 
a  Poisson  process. 


2.2  Product  moment  score  statistics 

An  alternative  to  the  serial  correlation  coefficient  is  obtained 

(Cox  and  Lewis,    1966,    pp.  166-7)    by  replacing  the  actual  interval 

values    x.     by  their  ranks     r.     or  exponential  scores      e(r.;n)     (Cox 
i  l  l 

and  Lewis,    1966,   pp.  26-27),    and  computing  a  first  order  product 
moment  statistic 


n-1 

R(n)  =    S     e(r.;n)X    e(r       ;n).  (2.9) 

1  i=1         l  l+l 


A  test  for  serial  independence  is  then  based  on  the  null  distribution  of 
R  (n)    under  a  permutation  hypothesis. 

An  advantage  to  using    R.(n)    is  that  it  controls  outliers  in  the 
series  of  events;    i.  e.  ,    missing  points.     Its  distribution  has  been  tabu- 


24 


lated  by  Lewis  and  Goodman    (196  9,    197  0)  for  ranks,    exponential  scores, 
scores  from  gamma  populations  with  parameter    k  =   -  1/2    and  scores 
from  Weibull  populations  with  shape  parameter    1/2     (1/2  Weibull). 

The  distributions  for  these  product  moment  rank  statistics  are 
similar  to  those  for  the  ordinary  serial  correlation  coefficient    p      from 
the  equivalent  parent  population;    positively  skewed  and  very  slow  to 
converge  to  the  asymptotic  normal  distribution.  The  distribution  of  the 
normalized  exponential  score  product-moment  statistic  does  in  fact 
converge  to  the  distribution  of    p      when  the  series  length  is    n»  10,  000; 
quantiles  are  shown  for    n  =  45  0    in  Table  5  and  should  be  compared 
with  the  exponential  values  in  Table  4. 

Table  5.      Normalized  exponential  score  product-moment  statistic;  n  =  450 

n*/  ess  r*s  ****  cw  rss  ^*y  /**» 

X0.  001        X0.  002        X0.  005        X0.  010        X0.  020        X0.  025        X0.  050        X0.  100 
-2.799         -2.623         -2.370        -2.159         -1.929        -1.848        -1.567         -1.240 


X0.  900        X0.  950        X0.  975        x0.  980        X0.  990        X0.  995        x0.  998        X0.  999 
1.270  1.656  1.994  2.096  2.395  2.669  3.009  3.245 


The  sampling  error  in  the  values  in  Table  5  is  approximately    0.  001. 
Note  that  the  idea  behind  these  tests  is    similar  to  that  for  the 
technique  of  random  shuffling  described  in  Perkel  et  al  (1967a)  and 
apparently  commonly  used  in  neurophysiological  work.      The  use  of  the 


25 


product -moment  statistic  is  more  sophisticated  and  does  not  require  a 
computer  to  do  the  shuffling. 


2.  3  Tests  of  serial  independence  based  on  the  periodogram 

The  first  serial  correlation  coefficient    p      can  be  used  as  a  test 
for  serial  independence  for  alternatives  other  than  the  first  order 
Markov  interval  process,    but  has  certain  drawbacks.     In  particular, 
for  the  more  common  alternatives  such  as  clustering  processes,    there 
is  no  imperative  for  looking  at  just  the  first  serial  correlation  and  tests 
combining  serial  correlations  of  several  orders  present  difficulties  be- 
cause of  the  correlation  between  these  statistics. 

Tests  based  on  the  periodogram  were  advocated  by  Bartlett  (1954) 
and  described  in  Cox  and  Lewis     (1966,  pp.   168-17CJ  and  Durbin  (1969). 

Denote  the  finite  Fourier  transform  of  the    n    intervals    x, ,  .  .  .  ,x      by 

1  n 


1  n  18tJp 

H    (u    )  =    ttz         2     x     e  ,  u)     =  2Trp/n, 

n      P  ,^^Z        8=1      8  P 


(27rn) 
and  the  periodogram  by 


I   (w   )=     |H    (w   )|2. 
n      p  n      p 


=  1,2,..,  [1/2  n]   =  I. 


The  test  is  essentially  for  a     "flat"     spectrum  using  the  distribution 
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theory  for  the    I   (w  )'s.      The  theory  states  that  the    I    (w   )'s    are 
n      p  n      p 

approximately  independent,    exponentially  distributed  random  variables, 
the  result  being  exact  for    x  ,  .  .  .  ,x      independent  and  identically  normally 
distributed. 

This  hypothesis  of  a    "flat"     spectrum  can  be  tested  using  a  modi- 
fied homogeneity  of  variance  statistic,    or  the  cumulated  periodogram 
values 

i  i 

U    =      S    I   (w  )  /    S     I    (w   )  p=  1,2,...,  1-1 

1        p=l    n      P        p=l     n      P 

can  be  tested  as  order  statistics  from  a  uniform  distribution.      Bartlett 

showed  that  the  distribution  theory  is  essentially  independent  of    k   , 

the  third  cumulant  of  the  intervals  between  events,    but  it  can  be  shown 

to  be  sensitive  to    k„,     the  fourth  cumulant. 

4 

Again,    results  have  been  obtained  by  synthetic  sampling  by 
Goodman  and  Lewis  (19*72).      Thedistribution  of  the  maximum  values 
of  the  periodogram  is  shown  in  Table  6,    the  values  for  normally  dis- 
tributed   x  's    being  exact,    those  for  exponential  and  half  Weibull  being 
taken  from  a  large  simulation.      The  deviations  from  the  normal  case 
for  exponential    x. 's    is  small;    for    1/2  Weibull    x  's,     which  have  a 
coefficient  of  kurtosis  of    84.720,    the  departure  is  dramatic. 

In  Figure  5,    for  exponential  variaties,    we  give  a  computer  plot 
of  the  distribution  of  the  modified  Kolmogorov-Smirnov  statistic 
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KSTWO(N)  =   s/n    C     =    sjh        max  |U    ,   -  i/l 


in  the  form  of  computer  printed  plots  of  sixteen  quantiles  as  a  function 
of    n,     the  length  of  the  series,      (n  =  N    in  the  plot  and  the  statistic  is 
called    KSCTWO(N)).      The  quantiles  are  a  little  smaller  than  for  the 
normal  case  in  which  the  upper     (asymptotic)    5  percent    and    1  percent 
points  are    1.  35  8    and    1.  62  8,     respectively.     Note  that  for    N  =  140    the 
statistic  is  based  on    t  -  69      U      's,      so  that  the  convergence  of  the 
quantiles  is  faster  than  it  appears  to  be  in  the  figure,    although  not  as 
fast  as  in  the  normal  case.      Of  course,   it  is  not  known  that  the  dis- 
tribution based  on  the  spectrum  for  non-normal  variates  converges  to 
the  distribution  for  normal  variates,    but  it  is  likely  if  the  first  four 
moments  exist. 
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3.     TESTS  FOR  POISSON  PROCESSES 

We  will  discuss  here  primarily  tests  for  the  Poisson  hypothesis 

against  stationary  alternatives,  trends  being  discussed  in  Section  5.      We 

assume  as  before  that    n    events  are  observed  at  times  t ,<  t_  <  .  .  .    <  t 

1       2  n 

in  the  fixed  interval    (0,  t   ].     Intervals  between  events  are  denoted  by 

x.,  x_,  .  .  .  ,  x   ,     and  it  is  convenient  to  denote  the  residual  interval  at  the 
12  n 

end  of  the  period  of  observation  by    x         =  trt  -  t   . 
"  y       n+1        0        n 

There  are  four  major  categories  of  tests. 


a)  Tests  based  on  the    t.  's,      conditional  upon    n  =  N(t    ),  which  is  a 

sufficient  statistic  for  the  nuisance  parameter      \,     the  rate  in  a  homo- 
geneous Poisson  process.      Since  the    t. 's    are    (conditionally)     the  order 
statistics  from  a  uniform  distribution,    the  empirical  count  function    N(t) 
(see  Figure  1)    is  proportional    to  the  empirical  cumulative  distribution 
function  for  the  uniformly  distributed  samples.      The  intervals    x.     are 
then  just  the  gap  statistics  for  the  sample. 

Tests  based  on  the  maximum  deviation  between  N(t)  and    t/t 
(Kolmogorov-Smirnov  statistics  and  modifications)    or  other  metrics 
(Anderson-Darling  statistic)    have  well  known  distribution  theories    (see, 
for  example,    Durbin,    1967),     but  are  sensitive  mainly  to  trend  departures 
from  the  homogeneous  Poisson  hypothesis.     In  fact,    Lewis    (1965)     showed 
for  the  special  case  of  gamma  renewal  alternatives,    that  the  test  is  not 
consistent.      Those  results  can  probably  be  extended  by  results  on 
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empirical  processes  for  renewal  processes  reported  in  Pyke  (1972). 

The  most  important  thing  to  note  here  is  that  while  the  distri- 
bution theory  of  these  tests  is  well  known  because  of  the  identity  with 
the  problem  of  testing, in  a  random   sample  of  size    n,     a  distribution 
F    (t)     against  an  alternative    F  (t),     the  distribution  theory  under 
alternatives  is  completely  different.     In  the  one  case    F     (t.)     are  order 
statistics  from  nonuniform  random  samples;    for  point  process  other 
than    (homogeneous  or  nonhomogeneous)    Poisson  process,    both  the 
uniformity  and  the  randomness     (independence)     disappear. 

Note  that  the  dispersion  test  for  Poisson  variates     (Cox  and 
Lewis,    1966,   p.  158)    is  equivalent  to  testing  the  uniformity  of  the    t  's 
with  a  chi-square  goodness  of  fit  statistic     (Lewis,    1965). 

b)  The  gap  statistics  for  the    x. 's,     i  =  1,  .  .  .  ,  (n+1),     which,    under 

the  Poisson  hypothesis,    are  a  random  number    n    of  independent  exponen- 
tial variates,    are  formed  from  the  order  statistics  of  the    x  's    as 


D    .  =  (x...    -  x..       )(n  +  2   -  i)  (x.  =  0,     i  =  1,  .  .  .  ,  n+1).    (3.  1) 

ni  (l)         (i-1)  0 


These  are  again  a  random  number  of  independent  exponentials 
(Cox  and  Lewis,    1966,    p.  26)    with  sum    t   ,     and  the  statistics 
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i 
t!  =      ED. 

1      j=i     nJ 


=     {x(1)  +  x        +  .  .  .    +  (n  +  2   -  i)X(i)  }/tQ      (i  *  1,  .  .  .  ,  n)         (3.  3) 


are,    conditional  upon    n,     order  statistics  from  a  uniform     (0,  1)    distri- 
bution   (see  Cox  and  Lewis,    1966,   p.  154    for  a  more  complete  derivation; 
the  conditioning  upon    n    and  the  fixed  total    t      is  very  subtle). 

A  great  deal  is  now  known  about  the  null  distributions  of  test 

statistics  based  on  the  gaps    D    .,     or  on  the  ordered    x. 's,     and  the 

6   K  ni  l     ' 

reader  is  referred  to  the  excellent  review  articles  by  Pyke    (1965,    1972). 
Also,    against  renewal  alternatives,    a  good  deal  is  known  about  asymptotic 
and,    in  some  cases,    small  sample  power.     Alternatives  are  generally 
specific  distributions  on  those  in  which  the  distribution  is  limited  by 
specifying  that  the  intervals  have,    say,    distributions  with  monotone  increas- 
ing or  decreasing  failure  rate    (Proschan  and  Pyke,    1967).     Again  see 
Pyke    (1972)    for  a  summary,    and  specifically  Jackson  (1967),    Bickel 
(1969),    and  Bickel  and  Doksum  (1969). 

Empirically  and  from  a  common  sense  point  of  view  it  is 
quite  clear  that  against  stationary  alternatives  tests  for  Poisson 
processes  based  on  the    t."s    are  more  useful  than  tests  based  on 
the    t.'s.      This  has  been  observed  in  using  the  SASE  program, 
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in  which  the  uniformity  of  the    t."s    is  tested  using  both  one-sided  and 
two-sided  Kolmogorov-Smirnov  statistics  and  the  Cramer-von  Mises 
statistic  (Cox  and  Lewis,    1966,   p.    147).  , 

More  information  about  the  power  of  these  procedures,    especially 
against  cluster  alternatives  would  be  valuable  (Lewis,    1969;  Vere-Jones, 
1970). 

c)  There  are  also  specific  tests  against  renewal  alternatives. 
The  Moran  test    (Cox  and  Lewis,    1966,   p.  161)    is  an  example. 

d)  A  useful  test  can  be  based  on  the  empirical  spectrum  of  counts, 
as  pointed  out  by  Bartlett    (1963).      Thus,    let  the  finite  Fourier -Stieltjes 
transform  of  the  sample  function    N(t)    be 


t 

H     (W)=  (irt   )"1/2       f°      eitU  c 
0 


£ 

^t=0 


1/2       n  n 

=  (irt J  {    S    cos(t.w)  +  i    2    sin(t.u)}  (3.5) 

0  j=l  j  j=l  j 


=  (irt  )"1/2{A    («)  +i  B     ((0)},  (3.6) 

0  0 


and  the  periodogram 


I     (u)  =    t«tftrl*{Af  (w)  +B2  (w)}. 
0  u  Z0  t0 
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We  will  refer  again  to  the  distribution  theory  of    I     (u>)     later. 

fc0 
Noting,     however,    that  for     wt     =  p2ir ,     I     (w)    is  approximately 

0  fc0 

exponentially  distributed  with  mean       X  tt  ,     and  that  for  any  two  such 

frequencies  the  correlation  between  the  periodogram  points  is     ~l/(l+n), 

we  can  see  that  the  spectral  test  for  independence    (Section  2)     can  be 

applied  here  to  test  for  a  Poisson  process. 

The  main  drawback  here  is  that  while  the  spectrum  of  intervals 
is  limited  to  approximately    n/2       periodogram  points,    the  count  spectrum 
is  not  limited  in  such  a  way.      Since  there  are  roughly     n/2      degrees  of 
freedom  available,    it  is  a  problem  as  to  which     n/2      points  of  the 
periodogram  with  frequency  of  the  form     u)t     =  2-rrp   to  use.      Using  more 
would  invalidate,   for  example,   use  of  standard  distribution  theory  of 
Kolmogorov-Smirnov  statistics  to  test  the  cumulated  periodogram. 

This  is  a  point  that  needs  considerably  more  work.      The  tests 
are  still  useful  in  an  informal  way,    particuarly  since  the  shape  of  the 
spectrum  can  suggest  physical  reasons  for  departures  from  a  Poisson 
process.     In  neurophysiology    (Pe  rkel  et  al,    1967),    the  tendency  has 
been  to  use  the  estimated  intensity  function    (Cox  and  Lewis,    1966,    p.  121) 
rather  than  the  spectrum  of  counts  to  assess  departures  from  the  Poisson 
hypothesis.      This  is  because  many  of  the  neurophysiological  processes 
causing  departures  are  more  simply  expressed  in  time  than  in  frequency. 
However,    the  distribution  theory  for  the  estimated  intensity  function  is 
difficult.     Cox  (1965)  has  discussed  the  distribution  theory  for  the  Poisson 
case. 
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4.      STATISTICAL  ANALYSIS 

OF  NON-RENEWAL  POINT  PROCESSES 

We  consider  here  several  non-  renewal  point  processes  which  are 
of  great  practical  importance.      For  the  most  part,    however,    the  struc- 
ture is  too  complicated  to  write  down  a  likelihood  function  so  that 
estimation  and  testing  is    ad  hoc.      Examples  of  such  analyses  are  cited. 

4.  1  Cluster  point  processes 

Cluster  point  processes  (branching  point  processes)  are 
important  because  they  arise  naturally  in  practice  and  have  in- 
teresting mathematical  properties  (Lewis,    1969;  Vere-Jones, 
1970). 

Briefly,    a  main  process     (usually  a  Poisson  process)     generates 
at  each  point  a  sequence  of  subsidiary  events.     In  Vere-Jones  (197  0)  the 
main  processes  are  initial  earthquake  shocks,    the  subsidiary  or  cluster 
events  are  aftershocks.      The  processes  of  subsidiary  evens  are  inde- 
pendent and  can,    in  general,    be  arbitrary  point  processes  which 
terminate  with  probability  one  after  a  finite  number  of  events  occur. 
Two  important  special  cases  arise: 

a)  The  subsidiary  events  are  generated  cumulatively  as  a  finite 

renewal  case.      This  is  known  as  the  Bartlett-Lewis  process. 
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b)  The  subsidiary  events  are  generated  additively,    each  one  of  the 

random  number  of  events  being  independently  displaced  from  its    (main) 
generating  event.      The  resultant  subsidiary  process  is  generally  non- 
stationary  and  the  complete  process  is  called  a  Neyman-Scott  cluster 
process. 

Computer  failure  patterns  generated  by  this  type  of  mechanism 
have  been  analyzed  by  Lewis  (1964),    and  earthquakes  by  Vere-Jones  (1970), 
Vere-Jones  and  Davies  (1966),    and  Shlien  and  Toksbz  0970).     A  fairly  good 
ad  hoc    analysis  can  be  given  for  the  Bartlett-Lewis  process  since  the 
marginal  distribution  of  intervals  is  known  (Lewis,    1964),    as  well  as  the 
spectrum  of  counts  and  in  some  cases  the  spectrum  of  intervals     (Gilles 
and  Lewis,    1967).      The  Neyman-Scott  process  does  not  yield  a  simple 
expression  for  the  interval  distribution  and  it  is  not  yet  known  if  the 
coefficient  of  variation  of  the  intervals  is  greater  than  one,    as  it  is 
for  the  Bartlett-Lewis  process. 

These  cluster  processes  are  overdispersed  relative  to  a  Poisson 
process  and  the  variance  time  curve  has  an  asymptotic  form  which  is 
independent  of  the  fine  structure  of  the  subsidiary  processes.     This  is 
a  help  in  analyzing  the  data;     basically  for  large  time  periods  the  counts 
of  events  N(t)  behave  as   though  the  subsidiary  events  were  concentrated 
at  the  main,    generating  event,     i.  e.  ,     like  a  bulk  Poisson  process. 
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Despite  these  points,    the  situation  with  regard  to  cluster  processes 
is  unsatisfactory;     questions  such  as  the  discrimination  of  Bartlett-Lewis 
and  Neyman-Scott  processes  arise  in  practice  and  are  not  solved    (see         ' 
the  discussion  in  Vere-Jones  (197  0)). 

When  the  cluster  process  has  clusters  with  only  one  event,    the  pro- 
cess is  equivalent  to  an  infinite  server  queue.     For  Poisson  main  events 
(parameter    \)     the  process  of  subsidiary  events  is,    in  equilibrium,    a 
Poisson  process  and  the  delay  distribution  cannot  be  determined. 

An  important  special  case  occurs  when  the  main  process  is  regu- 
lar and  each  event  is  independently  delayed.     Appointment  processes  are 
very  often  of  this  type  and  the  determination  of  the  delay  distribution  from 
observation  of  the  subsidiary  events     (arrivals)    has  been  considered  in 
detail  by  Govier  and  Lewis  (1967)    for  several  delay  distributions.      Since 
the  variance-time  curve  has  a  finite  limit,    they  have  called  these  pro- 
cesses with  controlled  variability. 

4.2  Superposed  processes 

Statistical  analysis  of  processes  which  are  superpositions  of 
point  processes  is  important,    particularly  in  neurophy Biological  contexts, 
and  usually  hinges  on  questions  of  estimating  the  number  of  contributing 
processes  or,   when  this  is  known,   identifying  the  structure  of  the  com- 
ponent processes. 
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These  are  difficult  questions  whose  solution  has  not  progressed 
much  beyond  the  basic  work  of  Cox  and  Smith  reported  in  Cox  and  Lewis 

(1966,    Chapter  8).      Without  specific  assumptions;    i.  e.  ,    that  the  com- 

i 

ponent  processes  are  renewal  processes  with,    say,    gamma  distributed 
intervals,    very  little  can  be  done  to  estimate  the  number  of  processes. 
The  problem  of  Walloe  et  al  (1969)  is  of  this  kind,    the  nature  of  the 
input  processes  and  the  neural  mechanism  being  well  known,    the  main 
question  being  how  many  inputs  impinged  on  the  neuron. 

Identifying  the  component  processes  is  again,    difficult,    without 
specific  assumptions.      Work  such  as  that  of  Ross  (197  0)  on  identifying 
the  interval  distribution  in  a  known  number    k    of  superposed  processes 
is  rather  technical;    note  that  the  variance  time  curve  and  spectra  of 
counts  are  additive;     e.  g.  ,    the  spectrum  of  counts  of  the  superposed 
process  is    k    times  the  spectrum  of  the  individual  processes.     Thus,  since 
the   spectrum  of  counts  of  a  renewal  process  is  simply  related  to  the  Laplace 
transform  of  the  probability  density  of  the  intervals,     f*(s), 


x 


.  f*(iw)  f*(-iw) 

«,<«>-(*»*  ft  +  ITf*^    +    ^7^7)  >  <4'» 


where     p.    is  the  mean  of    x,     the  spectrum  of  the  superposed  process 
g    (u))  =  kg    (u)    and  it  is  natural,    and  probably  efficient,    to  use  the 
estimated  spectrum  or  the  estimated  covariance  density  to  estimate 
f   (x).      Convergence  follows  from  results  of  Brillinger  (1972). 
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Nonhomogeneous  superpositions  also  occur     (Blumenthal,    Green- 
wood,   and  Herbach  (1971)),    but  as  in  the  homogeneous  case,    when  the 
number  of  contributing  processes  is  large,    the  resultant  process  over 
realistic  periods  of  observation  is  almost  indistinguishable  from  a 
(homogeneous  or  nonhomogeneous)    Poisson  process. 


4.  3  Doubly  stochastic  point  processes 

Let    A(t)     be  a  real  valued,    nondecreasing  stochastic  process. 
A  doubly  stochastic  point  process  is  a  generalization  of  the  nonhomo- 
geneous Poisson  process  in  which  the  integrated  intensity  function  is 
replaced  by    A(t).      Thus,    given  a  realization  of    A(t),     the  point  process 
is  a  nonhomogeneous  Poisson  process.     Generally    A(t)    is  differentiable 
and    X  (t)  =  A'(t)    might  be  a  stationary  stochastic  process,    a  determin- 
istic trend  or  a  combination  of  both.      The  process  was  introduced  by 
Cox  (1955);     see  Bartlett  (1966,    p.  325),    Cox  and  Lewis  (1966,   p.  179), 
and  Grandell  (1970)  for  details.     Gaver  (1963)  considered  the  case  where 
X  (t)     changes  level  at  random  times  and  called  the  process  a  random 
hazard  process. 

The  doubly  stochastic  mechanism  in  a  point  process  is  very 
realistic  and  probably  quite  common.      Thus    computer  failure  processes 
depend  to  some  extent  on  temperature,    humidity,    etc.      Unfortunately, 
analytic  properties  of  the  process  when    X  (t)    has  a  stochastic  element 
are  very  difficult  to  derive;     see  Lawrance  (1972).     If    X  (t)    is  a  stationary 
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-  2 

stochastic  process  with  mean     \,     variance    <*        and  autocorrelation 

A. 

function 


co  -lux 

p(t)  =  1/?tt     /        p(t)  e     W    dF(u),  (4.2) 

-  CO 


where    F(u)    is  the  integrated  spectrum  and    f(u))  =  F*(u)    when  it  exists, 
then  several  useful  results  can  be  obtained    (Cox  and  Lewis,    1966, 
p.  179-183).     In  particular,    the  covariance  density     v    (t)    of  the  doubly- 
stochastic  Poisson  process  is 


V+ft)  =  °l    p(t),  (4.3) 


and  the  spectrum  of  counts 


g    (U)  =    X/tt     +  2f(u)).  (4.4) 


Also,    in  general,    the  index  of  dispersion  is 


I(co)  =  1+2   &Z       f°°  p(u)  du,  (4.  5) 

X  0 


so  that  the  process  is  overdispersed  relative  to  the  Poisson  process. 

There  has  been  a  great  deal  of  interest  in  the  doubly  stochastic 
Poisson  process  in  physics,    optics,    and  engineering,    and  some  work 
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on  the  estimation  problems. 

a)  In  physics  photon  emission  is  a  well  understood  physical  process, 

and  it  is  known  that  the  process  modulating  the  Poisson  emission  is,    for 
instance,    an  Ornstein-Uhlenbeck  process.      The  parameters  of  this 
process  are  physical  parameters  which  it  is  of  interest  to  estimate, 
and  attempts  have  been  made    (e.  g.  ,   Pusey,    1971,   Jakeman,  Pyke  and 
Swain,    1971,  and  Koppel,    1971)  to  do  this  via  the  spectrum  and  the 
relationship  (4.  4).     Again,   laser  light  is  deliberately  focused  on  a  physical 
system  and  the  resultant  intensity  fluctuations  in  the  scattered  light  re- 
flects rates  of  molecular  motions  and  interactions  in  the  system. 

In  most  cases  the  number  of  photon  counts  in  these  experiments 
is  very  large  and  it  is  convenient,    and  very  often  necessary  to  cumulate 
counts.      There  are  problems  of  determining  the  best  sampling  interval, 
a  problem  which  is  often  made  simpler  by  detailed  knowledge  of  the 
modulating  process. 

Perhaps  because  of  the  large  amount  of  data  involved  and  the 
fact  that  the    \  (t)    process  is  fairly  well  understood,    physicists  rarely 
bother  about  details  such  as  the  efficiency  or  optimality  of  estimation 
procedures.      Moreover,    the  spectral  estimation  is  very  often  automated 
and  done  digitally. 
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b)  In  engineering,    two  areas  are  involved.      One  is  biomedical 

engineering  in  which,  for  example,    radioactive  substances  are  injected 

into  the  blood  stream  and  the  counts  of  the  radioactive  emission  are  used 

i 
to  estimate  a  decay  function  which  is  related  to  a  physical  process  of 

interest    (Snyder,    1971). 

A  second  area  is  optics,    where  modulated  light  is  used  to  trans- 
mit information    (Reiffen  and  Scherman,    1963;     Bar-Davidf    1969;  Karp 
and  Clark,    1970;     Clark  and  Hoversten,    1970).     Here    X  (t)    is  very  often 
a  signal  changing  levels,    a  possible  problem  being  the  discrimination 
of  several    X  .(t)'s    from  photon  counts  of  the  noise.     Very  often  a 
stationary    "noise"    element  is  also  present;     see  Bedard  (1966)  for  the 
physical  considerations. 

In  all  of  this  engineering  literature,    there  is  much  concern 
with  optimality,   in  some  sense,    and  procedures  are  based  on  likelihood 
ratios,    Bayesian  posterior  statistics  and  maximum  likelihood  detectors. 
It  is  a  difficult  literature  to  penetrate  and  my  overall  impression  is  that 
there  are  many  hidden  assumptions  involved,    one  being  that  samples 
are  large,   the  other  a  normality  assumption  which  may  be  quite  incorrect. 
Moreover,   most  explicit  results  are  for  very  simple  situations  such  as 
mixed  Poisson  processes  of  one  type  or  another,   or  rates  (Reiffen  and  Sherman, 
1963;  Snyder,    1972)  changing  at  known  times.     The  latter  problem  is  very  simple 
and  straightforward,    e specially  when  compared  to  the  case  of  unknown  change 
points  which  gives  a  true  doubly  stochastic  Poisson  process  and  an 
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inference  problem  similar  to  the  difficult  change  point  problem. 

c)  Grandell  (197  1,   1972)  has  considered  inference  in  doubly  stochastic 

processes  from  the  viewpoint  of  mathematical  statistics,    his  motivation 
apparently  being  primarily  problems  in  actuarial  work.     His  approach  is 
quite  different  from  that  in  the  engineering  literature,    cited  above,    and 
he  has  considered,    for  example,    optimum  estimates  of    E{\(t)}    based 
on  linear  combinations  of  the  data.      The  problems  are  related  to  the 
general  problem  of  curve  estimation  from  random  data  and,    while 
Grandell  (1972)  gives  specific  examples,    much  remains  to  be  done. 
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5.       TREND  ANALYSIS 

IN  NON-HOMOGENEOUS  POISSON  PROCESSES 

We  discuss  now  the  analysis  of  trends  in  non-homogeneous  Poisson 

processes,    extending  the  results  of  Cox  and  Lewis  (1966,    Ch.    3).     Initially! 

we  discuss  results  based  on  specific  parametric  models  for  the   rate  function 

X(t)    of  the  non-homogeneous  Poisson  process.     These  results  are  based  on 

the  fact  that  the  likelihood  function  for    n    observations  in  the  fixed  period 

(0,  t    ]    at  times    t<t^< <t      is 

'    0  12  n 

fco 

L(t. t   ,n;»)=  H     \(t  ;  9)    exp  {  -  /      \(u;e)du},  (5.1) 

1  n       ~~      1  =  1         i   —  u  — 

where  ,Q     denotes  the  vector  of  parameters  in  the  model.     Moreover,    given 
that    n    events  occur  in     (0,t    ],     the  times  to  events    t      are  the  order  statis- 
tics from  a  random  sample  from  the  probability  density  function 


f(t;e)=-M^) =  ^1SL,    0.t5t0,  (5'2) 

J00\(u;9)du  A(t0;§) 

and  the  conditional  likelihood  is 


n 


n 


L(t.,...tt   ;n;6)  =n»    II.  \(t.;9)/{A(tn;9)}    .  (5.3) 

Later  in  the  section  we  discuss  procedures  for  examining  the  adequacy 
of  the  model  for    K(t)    and  adequacy  of  the  non-homogeneous  Poisson  process 
model  itself. 
5.  1  Monotone  and  evolutionary  trends 

The  estimate  of    \(t)    when  there  is  no  trend  present  i.e.   X(t)  =  K,     is 
n/tQ. 
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In  Cox  and  Lewis  (1966,.  Chapter  3)  the  exponential  linear  trend 

Mt)  =  exp  {a+£  t}   =  Xexp{0  t}  (5.4) 

was  discussed  and  it  was  shown,    using  (5.  1),    that    (n,Zt  )    were  a  set  of 

*  i 

sufficient  statistics  for    («,0).     A  test  for   0=0    against   0^0,     where    a    is 
a  nuisance  parameter,    is  based  on  the  distribution  of    Zl./n,      given    n.     This 
is  an  optimum  (conditional)  test.     The  maximum  likelihood  estimate  of    /3 
was  also  given  implicitly,     but  no  small  sample  properties  of  the  estimate 
are  known. 

The  test  for    fi  =  0    was  applied  in  Section  1   to  analyze  a  sequence  of  ar< 
rivals     at  an  intensive  care  unit  (Table  1)  and  a  sectioned  analysis  indicated 
that  the  trend  was  not  monotone,    although  same  overall  increase  in  the    rate 
was  present. 

An    exponential     quadratic  trend  uses  the  model 

2  2 

\(t)  =exp{*+0t+  \t   }=  \exp{0  t+  yt   },  (5.5) 

for  which  we  get 

lnL(<r,0,<Y)  =  nln  X  +0Zt   +-Yi;t  (5.6) 

'<>  2 

-\/Q  exp  (/3u+yu    )du. 

The  maximum  likelihood  estimates  of    X,  o,/5     are  given  as  the  solution  of  the 

equations 

t 
t  =  n/  {/0°«p  £u+yu2)du},  (5.7) 

2?t  /  °u  exp($u  +\u2)du 

— - L  -      0 =0  (5   8) 

n  /o         ,■*■       +   2V 

J     exp(/3u+\u    )du 
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^.2  So    2  -         -     2,  , 

Zt?  /      u    exp  0  u  +  "Y  u   )  du 

—     "        J i    0,  (5.9) 

0  2 

/     exp(£u+^u    )du 

2  ' 

and  it  is  clear  that    (n,     Zt.,    Zt    )    are  a  set  of  sufficient  statistics  for 

(«t  0>   \). 

There  are  several  interesting  open  questions  here. 

What  are  the  small  sample  properties  of  the  estimators  and  what  are 
their  large  sample  properties?     Note  that  the  usual  theory  for  maximum 
likelihood  estimates  is  not  directly  applicable  because  of  the  random  sample 
size. 

What  effect  does  the  quadratic  term  (i.e.   "Y  /  0)    have  on  the  estimates 
of   0     in  the  model  (5.  4)  ? 

For  the  arrival  data  the  following  estimates  were  obtained. 

Linear  model::  X.  =    0.6342;}  j}  *  0.  0001427;  ^  , 

Quadratic  model:  \  *    0.4792;  p  =  0.  0009788:  \=~0.  4481x10 

Note  the  large  difference  in   /3     in  the  two  cases. 

Another  problem  of  Interest  is  to  test  for  the   quadratic  term  in  the 
trend.     It  is  clear  from  Table  1  that  either  a  quadratic  team  or  a  long  cycle 
is  necessary  to  adequately  model  the  annual  data  from  Section  1. 
A  formal  test  can  be  based  on  the  idea  that,   for  any  given    \,     n    and    Zt 

are  a  set  of  sufficient  statistics  for    a,/3.     Therefore  a  test  for    y    is  based 

2 

on  the  statistic    Zt     and  its  conditional  distribution,   given    n    and    Zt  .     This 

2 
distribution  is  difficult  to  obtain  analytically,     Zt.    being  the  square  of  the 

distance  to  the  sample  point,   which  is  constrained  to  lie  in  the    (n-1 )- dimen- 
sional hype rp lane  defined  by    Zt    =  c. 
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For  large  samples  this  distribution  can  be  obtained  from  the  fact 

2 
that    Zt ,/n    and    Zt    /n,      conditioned  on    n,     are  jointly  normally  distributed 

for  large    n,      and  the  following  exact  moment  results: 

i 

\i%  =    EfZtj/n}   =  tQ/2;         irj  =  v*x{  Z.tj/a}  =  tQ  (Ifcif   I 

u2  =    E{  Zt2/n}  =    tQ/3;         o-2  =  var  {  Zt2/n}  «  4t*  (45nf  * 

p=  corr{Bi,    Bi)  =    ±J  =  0.968. 
n  n  4 

2 
Thus,   using  normal  theory  results  we  test  with    Zt   /n    having  a  normal 

distribution  with  mean  and  deviation 

=  'o/3  +to(— l   "  12  ).  <5-10> 

n  2 

0-  =<r1   (Up2)*.  (5.11) 

For  the  complete  set  of  arrival  data  Zt  /n  =  954.  25,    giving   \i  =  1,  187,  783, 
and    0- =  6,  530,     Z  tf/n  =  1,  1  61,  565,    giving    (  Zt2/n-  \i)/<r  =  -4.  015    and    this 
is  clearly  a  highly  significant  rejection  of  the  hypothesis  that    "y  =  0. 

Boswell  (19  66)  and  Boswell  and  Brunk  (1969)  have  considered  tests 
of  the  hypothesis  that  X(t)  is  not  constant  but  is  non-decreasing.  Using  a 
likelihood  ratio  criterion,    conditional  on  the  fixed  value  of    N(tQ)sn,     he 


found  the  test  statistic 
n 

where 


(n,    t  T1  .  (5.12) 
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the  null  hypothesis  being  rejected  for  large  values  of  the   statistic. 

The  statistic  (5.  12)  is  not  simple  to  compute,    but  Boswell  gave  an  ' 
iterative  procedure  for  the  computation,    and  same  results  on  the  limiting 
distribution  of  the   statistic. 

It  would  be  of  interest  to  compare  the  power  of  this  test  against  the 
power  of  the  test  based  on     Zt   /n    for  the  exponential  linear  trend.     Some 
results  have  been  obtained  by  Mr.    Ian  White  of  the  University  of  Edin- 
burgh. 

5.  2    Cyclic  trends  of  fixed  frequency 

Cyclic  time  trends  (as  opposed  to  cycles  on  serial  number)  occur 
frequently  in  point  processes  but  were  treated  only  as  an  exercise  in  Cox  & 
Lewis  (1  966).     An  example  of  such  a  series  is  given  in  Forrest  (1950),   who 
was  investigating  thunderstorm  severity  in  Great  Britain  and  its  effect  on 
power  lines.     He  found  a  tendency  for  thunderstorms  to  occur  in  the  morning 
(time  of  day  effect)  and  of  course  a  very  strong  time  of  year  (seasonal)effect. 

The  arrival  data  discussed  in  Section  1  might  be  expected  to  have 
such  effects,   although  a  time  of  day  effect  is  by  no  means  evident  since  there 
is   only    about  one  arrival  every  day  and  a  half.     In  Figure  6  we  show  a 
collapsed  plot  of  the  numbers  of  arrivals  in  the  successive  hourly  periods 
of  all  days  of  observation.     The  plot  should  be  compared  to  Figure  1 .  8a  in 
Cox  &  Lewis  (1966)  for  the  arrivals  from  4  February  1963  to  18  March  1963. 
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Figure  6.      Arrival  of  patients  at  an  intensive  care  unit. 
Collapsed  hourly  plot  of  number  of  arrivals  to  investigate  "time-of-day"  effect, 
Second  period  of  observation:     19  March   1964  to  8  February  1968.      Total  number 
of  arrivals  (including  ties):     1216.      Observation  time    t     =  1420  days;  m  =  0.  856. 
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Figure  7.      Arrival  of  patients  at  an  intensive  care  unit. 
Collapsed  daily  plot  or  arrivals.      All  arrivals  including  ties,      Solid  line 
is  2nd  period,    19  March  1964  to  8  February  1968  (1216  arrivals  in  1420 
days),      Dashed  line  is  complete  record.   4  February   1963  to  8  February 
1968  (1467  in  1829  days).       .Z     (n.  -  n)Z/(n)  =  12.  93,      Note  that 

4.095~12-6- 
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The  plot  for  the  complete  period  (4February  1963-6  February  1968)  in 
Figure  6  is  similar  in  shape  to  the  plot  for  the  earlier  period  and  no  formal 
statistical  test  is  needed  to  conclude  that  there  is  a  real  time-of-day  effect;. 
As  a  model  for  the  rate   X(t)    in  a  non-homogeneous  Poisson  process 
with  fixed  frequency    u)    =  2tt  /T       we  take 

X(t)  =  exp    {ar+k     sin     u)Qt  +  k,    cos    U)t}  (5.  14) 

=    X  exp  {k    sin(u)t+0)},  (5.15) 

"7  0    1  1  /» 

where    k   =  <k     +  k     f  ,     0   -  tan"     (  —  )  ,   X  =  exp(a)xl   (k), 
s         c  k  0 

8 


and    In(k)    is  a  modified  Bessel  function  of  the  first  kind  of  zero  order.     This 

reparameterization  is  used  because 

T 

/0°    exp    {k  sin  (  %p  +  6  )  }  dt  =    IQ(k)  ,  (5.  1  6) 

so  that  if    t   ,     the  total  period  of  observation  is  a  multiple  p   of   T     (here  one 

day),   we  have 

t    - 

A(tJ=    /  °   X(u)du  =  Xt     =  XpT  .  (5.17) 

v  0  0  0  0 

In  effect  we  have   separated  out  multiplicatively  a  linear  growth  from  a 
cyclic  effect  which  takes  place  only  within  periods  of  length     T  . 
The  reason  for  using  the  form  5.  14   rather  than,    say, 

X(t)  =  a  +  k    sin  (u  t+9)  (5.  18) 

is.  that  X(t)    must  be  positive,    and  to  achieve  this  with  (5.  18)  requires  k<a. 


51 


For     k  «  a    the  two  models  are  essentially  the  same.     In 
addition,    the  rate    (5.  14)     gives  simple  results  based  on  sufficient 
statistics,    this  being  closely  connected  with  the  fact  that  the  density 
(5.  2)     belongs  to  the  exponential  family  of  density  functions. 

An  additional  reason  for  using  (5.  14)  is  that  such  nonsymmetrical 
cycles  are  probably  more  natural  with  series  of  events,   possibly  because 
of  the  positivity  of  the  rate  function.     Professor  J,  W.  Tukey  has  suggested 
a  model  which  is  (5,  18)  squared;  this  could  be  useful  with  data  in  which  there 
are  a  large  number  of  arrivals  in  each  period,    so  that  regression  techniques 
using  the  square  root  of  the  numbers  of  arrivals  in  fixed  intervals  can  be 
used  (Cox  and  Lewis,    1966,    Chapt.    3).      This  model  has  also  been  used  by 
Fisher  (1953). 

In  Lewis  (1970)  the  following  results  were  given  for  the  non-homo- 
geneous Poisson  process  with  the  rate    X  (t)  given  in  (5.  14). 

Using  (5.  3)  we  find  that  the  observations  enter  into  the  likelihood 

function  only  through 

n  n 

n,  A     (coj  =    2    cos  unt.,  B     (uj  s    S      sin  cot.,        (5.19) 

t0     °        i=l  °  1  t0     °         1*1  °  x 

and  these  are  a   set  of  sufficient  statistics  for  the  parameters    cr,    k 

s 

and     k      in    (5.  14).      Note  that    A      (co  )     and    B     (oj  )     are  the  components 

c  'o     °  lo     ° 

of  the  periodogram     (3.  6). 

Maximum  likelihood  estimates  of    X  ,  9,     k    are 
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A     <«0) 

9=  tan"  (  b-k)}  <5'20) 

V   ° 


^    =  n/tQ,  (5.21) 


and    k    is  the  solution  of  the  equation 


2  2  Mk) 

(l/n){A      (co   )  +  B     k))rl^=    ip(k)f  (5.22) 

fc0     °  'o     °  I0(k) 


where      I  (k)     is  the  derivative  of    In(k),    and     i{j(k)     increases  monotonically 
from     0    to    1    as    k    goes  from     0    to    co.      This  later  fact  allows  one  to 
use  the  Neyman-Pearson  lemma  to  show  that,    conditional  on  the  observed 
value    n    of    N(t   ),     the  most  powerful  test  of    k  =  0    against    k?i  0    is 
based  on  the  statistic 

A2  (to  )  B2  (co  ) 

•    Vn)  f-r—  +  ~T—  }=  <-V"»  \  *i>-  (5-23) 

0  0  0 

Since,    for    k  =   0,     21     (to   )    has  asymptotically    (Cox  and  Lewis, 
1966,    Chapter  5)    an  exponential  distribution,    the  test  for    k  =  0    reduces 
to  accepting  the  hypothesis  at,    say,      a     5%    level  if 


21     (co  )  =    fc/ir)X(t   /n)=    £>/tt)X   m  .  (5.  24) 

0 


2  .,  2 

The  factor     6    arises  because    prob  \\     =  6}  =  .  95,     where    x       is  a 
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random  variable  with  a  chi-square  distribution  of  two  degrees  of  freedom. 
This  test  applied  in  Lewis  (197  0)  to  the  first  section  of  the  arrival 
data  gave  a  strong  indication  of  the  presence  of  the  cycle  at  a  period  of 
one  day. 

For  the  complete  record  of  arrivals  at  an  intensive  care  unit  dis- 
cussed in  Section  1  we  get  for  the  periodogram  at     p  =  1829     or     go     =  2tt, 
which  is  the  day  frequency,    since     t     =  1829     and    go     =  2-rrp/t 


21     (go  )  =  27.  094    and     6m/TT    =  1.  52. 


This  is,    not  surprisingly,    highly  significant.     For  the  first    409 
arrivals  the   results  are  (Lewis,     1970) 


21     (go  )  =  5.  0    and    6m/-rr  =  1.  10. 
0 


Thus  the  periodogram  component  is  increasing  roughly  in  proportion  to 
n,     as  it  should  for  a  true  cyclic  component. 

When  Lewis  (197  0)  was  written,    the  connection  of  this  model  for 
a  cycle  of  fixed  frequency  with  tests  for  directionality  on  the  circle  was 
not  realized.     In  fact,    the  conditional  test    (5.  24)    is  equivalent  to  testing 
that    n    observations  on  a  circle    (here  a    24-hour  clock)    have  the  von  Mises 
or  circular  normal  distribution    (Gumbel,    Greenwood     and  Durand ,    1953) 
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e*p{ksin(i+9)}  ,o<(<2f|.  (5.25) 

Ztr  I0(k) 


For    k  =  0,     this  is  a  uniform  distribution;     otherwise,    it  has  a  mode  at 
9,     the  vector  in  the  direction    9    being  called  the  modal  vector.     Green- 
wood and  Durand  (1955)  obtained  the  distribution  of  the  square  of  the 
resultant 


R2  =  (   S    cos  I.)2  +    (  S    sin  i.)2  (5.  26) 

i=l  i  i=l  i 


of    n     such  vectors  when    k  =  0,     generalizing  earlier  work  of  Pearson 
on  the  problem  of  random  walks  on  the  circle.      The  formal  analogy  of 
(5.  26)    with  the  periodogram     (3.  5)     should  be  clear.      Watson  and 
Williams     (1956)     generalized  these  distributional  results,    using  results 
for  sufficient  statistics,    and  found  the  conditional  p.  d.  f.      of  the  quantity 
on  the  left  of  equation  (5.  22),    which  we  denote  by    r,     as 


co 


f(n'r)=  rIQ(kr)    J         {JQ(u)}     JQ(ru)du/ {lQ(k) }   ,  (5.27) 


where    J   ( •  )    is  the  ordinary  Bessel  function  of  zero  order. 

It  is  not  all  apparent  that    (5.  27)    is  more  useful  for  computation 
than  the  generating  function,    which  is  given,    for  example,    for    k  =  0 
in  Cox  and  Lewis     (1966,    Chapter  5)    in  the  discussion  of  the  distribution 
of  the  periodogram  at  one  point      cj   ,     where    co      is  of  the  form 
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t  oo     =  2-n-p,     for    p    integer.     However,    Stephens     (1969a,    1969b)    has 

used    (5.  27)    to  tabulate,    among  other  things,    the  power  of  the  test     that 

k  =  0  against  alternatives   k  >  0.     The  most  complete  discussion  of  these  tests 

for  the  von-Mises  distribution  is  given  in  Watson  and  Williams  (1956).     The 

function   i|;(k)  is  tabulated  by  Gumbel,    Greenwood,    and  Durand  (1953).  Stephens 

(1969b)  has  also  discussed  tests  for    9  =  9    ,     and  joint  tests  for    k    and    9. 

It  is  clear  that  problems  involving  cycles  at  two  or  three  fixed 
frequences,      e.  g.  ,     to   ,    to  ,    go   ,     will  arise  in  analyzing  series  of  events. 
For  example,    in  the  problem  of  analyzing  the  arrivals  at  an  intensive 
care  unit,    a  time  of  week  effect  and  a  time  of  year  effect  are  distinct 
possibilities.      Surprisingly  enough,    there  does  not  seem  to  be  a  strong 
time  of  week  effect  in  the  data,    but  there  is  not  space  here  to  go  into 
this. 

Formal  tests  for  more  than  one  cycle  follow  from  results  in 

Cox  and  Lewis     (1966,    Chapter  5)     to  the  effect  that  at  two  different 

frequencies    to      and    go  ,     both  of  which  when  multiplied  by    t         are    2tt 

times  an  integer,    the  correlation  between  the  periodogram  values 

I     (to  )     and    I     (to  )       is    approximately     (1  +  n)      .      Thus,    a  test  for  a 

0  0 

time  of  week  effect    (to.),     based  on  the  conditional  distribution  of    I(to  ), 

given  the  values  of  N(t    )  =  n    and    I     (w_),     reduces  in  effect  to  an  inde- 

0  tQ 

pendent  test  at    to      based  on    (5.  23). 

For  the  arrival  data,  considered  over  a  period  of    1456    days, 

we  get  for    p  =  208    (to     *  2ir/7),     or  a  period  of  one  week,    the  value 

P 
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AZ  (co  )           BZ  (co   ) 
2              fcO     P                fcO     P 
21,   («L>-     7       ("— +     -T >     =    0-467. 

'o    p  fco  fco 

Thus  there  is  no  formal  indication  of  a  time  =  of=week  effect. 


5.  3  The  spectrum  of  counts  and  cycles  of  unknown  frequency 

In  the  previous  subsection,    we  showed  the  connection  between 
tests  for  a  cyclic  component  at  a  known  frequency  in  a  nonhomogeneous 
Poisson  process  and  the  periodogram,    the  periodogram  being  the  basis 
for  estimation  of  the    (second  order)     spectrum  of  counts    g    (co). 
Clearly,    one  might  want  to  look  for  cycles  at  unknown  frequency  and 
this  will,    intuitively,    be  based  on  the  spectrum  of  counts.     More  pre- 
cisely,   it  will  be  based  on  the  periodogram 


I     (co)  =  (1/irt    )    {AZ  (co)  +  B^  (co)}.  (5.28) 

0  0  0 


The  analogous  problem  in  ordinary  time-series  analysis  is  the 
classical  problem  of  hidden  periodicities,    discussed  at  length  by  Hannan 
(1970,    p.  463).      This  problem  has  not  been  tackled  in  point  processes, 
and  is  more  difficult  for  two  reasons.      First,    the  periodogram  points 
are  not  quite  uncorrelated,    as  we  saw  in  the  previous  section,    and  also 
the   spectrum  is,    in  theory,    not  limited  in  the  frequency  domain.      (There 
will  always  be  some  band  limiting  due  to  jitter;     see  Lewis,    1970.  ) 
Thus,    it  is  a  problem  how  to  use  the  distribution  theory  for  the  spectrum 
given  in  Bartlett  (1963)  and  Cox  and  Lewis  (1966,    Chapter  5)    and  how  to 
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pick  relevant  parts  of  the  spectrum  from  which  to  estimate  the  unknown 

frequency. 

As  an  example,    we    give  in  Figure  8  the  spectrum  of  the  arrival 

data  of  Section  1.      The  peak  for  the  day  effect  at    p  =  1829  (co     =  2tt)    is 

P 

evident,    but  it  would  be  difficult  to  read  into  this  spectrum  anything  but 
a  conclusion  that  it  is  a  nonhomogeneous  Poisson  process  with  a  single 
fixed  cycle    (day-effect).      There  is  possibly  a  harmonic  at    p  =  3658 
and  some  inhibition  at  low  frequencies,    but  all  the  points    (except  for 
p  =  1829)    are  within    1  percent    confidence  bands  for  individual  values, 
and  would  be  within  bands  for  the  maximum  of  the  periodogram  values 
(Bartlett,    1963). 

The  spectrum  of  counts  is  generally  a  useful  tool,    and  even  more 
particularly  so  for  non-Poissonprocesses,  but  has    found  very  limited 
acceptance  amongst  applied  workers.      This  is  partly  due  to  confusion 
with  the  spectrum  of  intervals,    with  which  neurophysiologists,  for  example, 
are  much  more  familiar. 

Lewis  (19*7  0)  has  given  a  heuristic  justification  for  the  spectrum 
of  counts  and  discussed  computation  and  smoothing.  Bartlett  (1967)  has 
discussed  finer  points  in  this  type  of  analysis. 

Finally,    Brillinger  (1972)  has  given  a  general  theory  for  spectral 
analysis  of  point  processes,    and  has  defined  higher  order  spectra 
(cumulant  spectra  of  order    k)    for  point  processes.      These  are  gener- 
alizations of  the  cumulant  spectra  of  order    k    of  continuous  time  series. 
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Figure  8.      Arrival  of  patients  at  an  intensive  care  unit, 
Spectrum  of  counts  for  the  complete  record:     4  February  1963  to  8  Febraury  1968, 
n  «=   1458  arrivals  in  a  period  of    t       =  18  29  days.      60  point  uniform  smoothing. 
Bands  are  0.  95  and  0.  99  confidence  regions  for  individual  values,     m  «  0.  7972. 
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Figure  9.      Arrival  of  patients  at  an  intensive  care  unit. 
Averages  of  successive  groups  of  20  inter-arrival  times  after  detrending  with 

ft  A  A  A  A  A~AA' 

t    =  \     exp  {  a  +   0  +  k  sin  ( 2  it  u  +  6  }  du,   where    a,    /3  ,    k,    6    are  maximum  likeli- 

^0 
hood    estimates,      n  «=   1458  arrivals  in  1829  days  from  4  February  1963  to 

8  February   1968,     Overall  average  after  detrending:     1.  0088.      Homogeneity 

of  variance  statistic  =  309.  07  (x2,,:   mean  7  1,   a  «=   11.  92), 
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Brillinger's  work  is  based  on  a  general  spectral  representation  for    N(t), 
the  counting  process  of  the  point  process,    and  provides  an  asymptotic 
theory  for  the  spectral  estimates  of  order    2     (periodogram)    considered 
in  this  paper.     Brillinger's  paper  is  far  too  extensive  to  do  justice  to 
here.      The  problem  of  whether  cumulant  spectra  of  order  higher  than    2 
will  be  useful  in  practice  remains  open  and  should  be  explored. 
Brillinger's  spectral  decomposition  may  provide  an  answer  to  whether 
spectral  analysis  can  be  useful  as  a  representation  in  real  systems;  i.e., 


neurons. 


5.  4  Mixed  models 

The  arrival  data  considered  in  Section  1  has  been  shown  to  have 
a  long  term  trend  which  can  be  represented  by  the  exponential  quadratic 
function    (5.  5)     and  a  strong  day  cycle.     A  combined  model  for  this  data 
could  be  a  nonhomogeneous  Poisson  process  with  rate    X  (t)     of  the  form 


2  , 

X.  (t)  =  expfa  +  pt  +  vt     +  k  sin(co   t  +  9)}.  (5.29) 


t 
This  is  difficult  to  handle  because     .A(t)  =   /    X  (u)du    is  not  tract- 

0 

able.     However,    in  the  context  of  the  arrival  data  the  evolutionary  trend 
changes  little  within  a  cycle  and  an  approximation  in  which  the  linear 
and  quadratic  terms  are  assumed  constant  over  the  period  of  each  cycle 
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yields  a  tractable  model.     The  estimates  of  k    and  6     are  then  simply  (5.  20) 
and  (5.  22),      y    and    p    are  estimated  from     (5.  8)     and    (5.  9)     and    X.     is 
given  by    (5.  7)     with    In(k)    multiplying  the  term  in  the  denominator. 
For  the  arrival  data,  we  obtain  the  following  values: 


\  =     0.47  92; 

k  =     0.  654; 

9  =    3.  519; 

p  =     0.  0009788; 

v  =   -0.  4481  X   10"    . 

Clearly,    mixed  models   such  as     (5.  29)     will  be  needed  in  many 
cases  for  the  analysis  of  real  data. 


5.  5  Residual  analysis  of  point  processes 

Many  observed  point  processes  have  rate  functions  which  are 
clearly  not  easily  representable  by  simple,    tractable  models  such  as 
those  discussed  in  the  previous  section.     For  example,    it  is  common 
in  observing  arrivals  of  jobs  in  a  computing  center  to  find    (roughly) 
the  rate  increasing  sharply  throughout  the  morning,    dipping  around 
lunch  time,    picking  up  slightly  in  the  afternoon,    and  then  dipping  off 
sharply  at  night.      The  situation  is  fairly  simple  in  this  example,    since 
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many  days  observations  are  available,    showing  similar  variations, 
and  careful  pooling  gives  a  smoothed  tabulated  estimate  of  the  rate. 
This  rate  is  needed  to  model  scheduling  algorithms  for  the  computer. 

In  general,    the  problem  of  smoothing  the  point  process  to  obtain 
an  estimate  of    X  (t)    is,    in  my  opinion,    open,    despite  the  work  on  doubly 
stochastic  Poisson  processes  cited  in  Section  4,    which  should  be  applic- 
able.    It  is  generally  connected  with  the  statistical  problem  of  curve 
smoothing,    and  will  not  be  considered  further. 

Both  in  the  case  where    X  (t)    is  obtained  by  smoothing  or  from  a 
specific  model,    say    exp  (or  +  (3t),     it  will  often  be  necessary  to  test  the 
adequacy  of  both  the  trend  model  and/ or  the  assumption  of  an  underlying 
nonhomogeneous  Poisson  process. 

For  example,    an  intensive  care  unit  is  a  service  facility,    and  to 
provide  adequate  service,    one  mast  know  something  about  the  arrival 
process;    it  might  be  very  regular  if  all  arrivals  are  from  surgical 
operations,    or  exhibit  clustering  if  all  arrivals  are  from  the  emergency 
room.      The  form  of  the  underlying  process  will  also  affect  the  test  for 
cyclic  trend  given  in    (5.  23),     the  attempt  to  estimate  the  spectrum  and 
filter  out  the  spike  usually  being  a  self  defeating  process    (Bartlett,    1967; 
Lewis,    1970). 

Both  formal  and  infoTmal  methods,    usually  graphical,    are  re- 
quired for  this  question,    analagous  to  the  similar  problem  in  regression 
analysis  which  is  usually  approached  by  examining  the  residuals  after 
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estimating  parameters  for  the  mean  value  in  the  model. 

While  keeping  in  mind  that  there  are  many  ways  to  approach 
this  in  different  situations,    the  following  proposal  seems  to  offer  a 
systematic  approach  to  examining  such  questions     (Lewis,    197  0).      It  is 
well  known  that  on  the  transformed  time  scale 


t   =    J"       \  (u)du,  (5.  30) 

0 


a  nonhomogeneous  Poisson  process  with  rate  function    X  (t)     becomes  a 
homogeneous  Poisson  process  of  rate    1.      Thus,    we  proposed  to  examine 
the  residual  process     (t.)    where 


t.  =    /  *  X(u)du,  (5-31) 


and    X  (t)     is   some  estimate  of    X  (t). 

There  are  a  number  of  ways  of  looking  at  this  transform. 

1)  If    X  (t)    is  very  smooth  over  periods  compared  to  the  mean  time 

between  intervals,    we  can  write 


X.=  t.   -  t.    .    ~   X(t.    ,)*,.  (5.32) 

l         i  l-l  l-l     l 


Then  a  logarithmic  transformation  reduces  this  to  standard  regression 
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models,    and  the  exponential  trend  functions  discussed  above  are  appealing 
since  we  have    (approximately)     standard  linear  regression  methods. 
These  are  considered  in  Cox  and  Lewis     (1966,    Chapter  3). 

Of  course,  in  data  such  as  the  arrival  data  of  Section  1,  this 
approximation  is  not  possible,  and  the  integral  transformation  must 
be  performed,    perhaps  numerically. 

The  real  problem  comes,  of  course,  when  X.  (t)  cannot  be 
assumed  to  be  smooth  over  suitably  short  intervals  and  no  simple 
parametric  form  for    \  (t)    is  available. 

2)  If  the  analysis  is  based  on  the  conditional  distribution  of  the 

t. 's,     given    N(t   )  =  n,     then  the  transformation    (5.  31)    is  seen  from   (5.  2) 
to  be  just  proportional  to  the  probability  integral  transform  for  the  density 
X-  (t)/A{t)     with  estimated  parameters. 

3)  There  is  also  a  possibility  of  looking  at    (5.  31)     as  a  filter  in 

the  spectral  domain. 

Formal  methods  for  analyzing  the    t.     are  difficult  and  will  be 
discussed  elsewhere.      This  is   simplest  for    X  (t;9)    when  sufficient 
statistics  for    9    are  available.      The  probability  structure  of  the    t. 
process,    given  the  sufficient  statistics,    is  independent  of    9_.      Of 
course,    in  the  simplest  case  of    X  (t)  =  X  ,     with    X   *  n/t0»     it:  is  wel1 
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known  that  the    X.'s     are  the  gap  statistics  for  independent  uniform 
random  variables  and  are  asymptotically  exponentially  distributed 
and  independent    (Pyke,    197?.).      For  small  samples,    they  will  gener- 
ally be  correlated  random  variables. 

Informal,    graphical  methods  are  useful  and  are  discussed  in 
the  context  of  the  arrival  data  given  in  Section  1.      The  first  section 
of  this  data  was  transformed,    using    (5.  31)    and  an  exponential  linear 
trend,    and  discussed  in  Lewis  (197  0).     No  test  for  Poisson  processes 
performed  in  the     SASE  IV    program  gave  significant  indication  of 
departure.      The  spectrum  of  the  transformed  data  is  shown  in  Lewis 
(1970). 

Clearly  the  (exponential)  linear  trend  plus  day  cycle  is  not  adequate 
for  the  complete  set  of  data,    as  we   saw  in  Section  1.      We  examine  its  adequacy 
further  here.     Figure  9  shows  a  plot  of  the  average  of  successive  groups  of 
20  X  's  after  detrending  with  the  (exponential)  linear  trend  plus  day  cycle.    The 
model  for   \(t)  is  clearly  inadequate.     Figure   10  shows  the  spectrum  of 
counts  after  the  linear  transformation.      There  is  no  spike  corresponding 
to  the  one  day  period  and  no  indication  of  departure  from  a  Poisson  process. 
Similarly,    the  tests  discussed  in  Section  3  for  Poisson  processes  do  not  give 
any  great  indication  of  departures.       However,    the  estimated  asymptotic  slope 
of  the  variance -time  curve  is  approximately  8,    and  since  this  quantity  is  re- 
lated to  the  rate  of  change  of  the  low  frequency  end  of  the  spectrum,    there  is 
again  an  indication  of  inadequacy  of  the  form  of  X.(t)  or  the  hypothesis  of  a 
nonhomogeneous  Poisson  process. 
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A  quadratic  term  was  added  to  the  rate    X  (t),     giving  the  trans- 
formation 


t 

1  ,A  A  A     2         A  A.  , 

t  .  =  /       exp   \a  +  (3u  +  ^u     +  k  sin(2Tr  u  +  9)  }du,  (5.  33) 

1         0 


where  the  values  of  the  estimated  parameters  are  given  in  Section  5.  4. 
The  plot  corresponding  to  Figure  9  is  shown  for  the    (exponential) 
quadratic  detrending  in  Figure  11.      This  and  other  results  give  a  much 
clearer  picture  of  a  possible  cycle  or  quasi-cycle  in  the  data  with  period 
of  about  a  year  and  a  half.      The  example  is  discussed  in  more  detail  in 
Lewis  (1971),    where  local  smoothing  is  used  to  estimate    \  (t)    for  this 
quasi-cycle  and  test  the  nonhomogeneous  Poisson  hypothesis. 
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Figure   10,      Arrival  of  patients  at  an  intensive  care  unit. 
Spectrum  of  counts  for  complete  record  after  detrending  with 


f1 
t    =   \    exp   {a 

Jo 


+  /3u+  ksin(2iTU+0)}du,   where    a,   /3,   k,  G    are  maximum 


likelihood  estimates,      n  =   1458  arrivals  in   18?9  days  from  4  February   1963  to 
8  February  1968.      Bands  are  approximately     0.95  and    0.99  confidence  regions 
for  individual  values. 
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Figure   11.      Arrival  of  patients  at  an  intensive  care  unit. 
Averages  of  successive  groups  of  20  inter-arrival  times  after  detrending  with 

exp   {a  +  y3u+  y  u  +  k  sin  (2tt  u  +  6)    }  du,   where  a,    fi  ,    y,   k,    0     are  maxi- 
mum likelihood  estimates,     n  =   1458  arrivals  in    t     =   1829  days  from  4  February 
1963  to  8  February  1968.     Overall  average  after  detrending: 
Homogeneity  of  variance  statistic:  (  X2  , :  mean  7  1,   tr  -  1 1.  92) . 
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6.  CONCLUSIONS 

Much  further  work  is  needed  in  the  statistical  analysis  of  series 
of  events.      Particular  areas  that  come  to  mind  are  the  analogues  of 
the  procedures  in  Section  5  for  modulated  renewal  processes.      Likeli- 
hood functions  can  be  written  down,    but  results  are  fairly  difficult  to 
obtain  from  the  expressions  (see  Cox  (197  2)  for  details).      Both  for  these 
processes  and  the  nonhomogeneous  Poisson  process  there  is  a  great 
deal  of  work  to  be  done  to  justified  estimation  and  testing  procedures  based 
on  random  sample  sizes.      Some  of  these  justifications  have  been  given  by 
Brown  (1972). 

Finally,   it  is  hoped  that  new  methods  will  be  developed  which  will 
allow  some  of  the  nonrenewal  point  processes  to  be  analyzed  in  a  more 
systematic  way. 
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tistics involved,  and  to  testing  the  functional  form  of  a  trend  in  a 
nonhomogeneous  Poisson  process,  as  well  as  the  point  process  model  itself 
A  survey  of  work  in  special  processes  such  as  cluster  processes  and 
doubly  stochastic  Poisson  processes  is  also  given. 
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